home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / gecp_r < prev    next >
Encoding:
Text File  |  1994-12-20  |  1.7 KB  |  72 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Gaussian elimination with complete pivoting.
  4.  
  5. // Syntax:      GECPL = gecp ( A )
  6.  
  7. // Description:
  8.  
  9. //      Gecp computes the factorization P*A*Q = L*U, where L is unit
  10. //      lower triangular, U is upper triangular, and P and Q are
  11. //      permutation matrices.  RHO is the growth factor. 
  12.  
  13. //    This file is a translation of gecp.m from version 2.0 of
  14. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  15. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  16.  
  17. //-------------------------------------------------------------------//
  18.  
  19. gecp = function ( A )
  20. {
  21.   local (A)
  22.  
  23.   n = A.nr; n = A.nc;
  24.   pp = 1:n; 
  25.   qq = 1:n;
  26.  
  27.   rho = 0;
  28.   maxA = norm(A[:], "i");
  29.  
  30.   for (k in 1:n-1)
  31.   {
  32.     // Find largest element in remaining square submatrix.
  33.     // Note: when tie for max, no guarantee which element is chosen.
  34.  
  35.     colmaxima = max( abs(A[k:n; k:n]) );
  36.     rowindices = maxi( abs(A[k:n; k:n]) );
  37.     biggest = max (colmaxima);
  38.     colindex = maxi (colmaxima);
  39.     row = rowindices[colindex] + k - 1;
  40.     col = colindex + k - 1;
  41.  
  42.     // Permute largest element into pivot position.
  43.  
  44.     A[ [k, row]; ] = A[ [row, k]; ];
  45.     A[; [k, col] ] = A[; [col, k] ];
  46.     pp[ k, row ] = pp[ row, k ]; 
  47.     qq[ k, col ] = qq[ col, k ];
  48.  
  49.     if (A[k;k] == 0)
  50.     {
  51.       break
  52.     }
  53.  
  54.     A[k+1:n;k] = A[k+1:n;k]/A[k;k];          // Multipliers.
  55.  
  56.     // Elimination
  57.     i = k+1:n;
  58.     A[i;i] = A[i;i] - A[i;k] * A[k;i];
  59.     rho = max( rho, max(max(abs(A[i;i]))));
  60.   }
  61.  
  62.  
  63.   L = tril(A,-1) + eye(n,n);
  64.   U = triu(A);
  65.  
  66.   P = eye(n,n); P = P[pp;];
  67.   Q = eye(n,n); Q = Q[;qq];
  68.   rho = rho/maxA;
  69.  
  70.   return << L=L ; U=U ; P=P ; Q=Q; rho=rho >>;
  71. };
  72.